﻿// mod09ga_prepare.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//使用C盘 OSGeo4W的GDAL库快速开发GDAL程序。
//处理云遥服务器上 ERA5 数据，将U、V风数据提取出来，然后风速，然后生成3dtileset和glb数据。
//文件名样例 2025042400_0p25_pressurelevels.grib
//数据105波段，每5个波段为一个气压层，每层的第4、5波段为u、v风。气压层顺序1,2,3... 为低气压（高空）到高气压（地表面）顺序。
//YY的ERA5数据左上角坐标为(-180.0 , 90.0) 而欧空局下载的是 (0.0, 90.0) ，本程序是针对YY的，所以注意这个问题。

#include <iostream>
#include "gdal_priv.h"
#include <string>
#include <vector>
#include <stdexcept>
#include <sstream>
#include <iostream>
#include <fstream>

#include "yy_write_gltf.h"
#include "blh2ecef.h"
#include "ecef_tiler.h"
#include <direct.h> //only for windows.
#include "ajson5.h"

using namespace std;

struct UvLayer {
    vector<float> _u;
    vector<float> _v;
    vector<float> _l;
};

//风速对应颜色
struct MeshRange {
    float valmin = 0;
    float valmax = 0;
    float r = 0.0;
    float g = 0.0;
    float b = 0.0;
    float a = 0.0;
};

vector<MeshRange> g_MeshRanges = {
    { 10.0,17.5, 0.33, 0.63, 0.70 , 0.2 }, //两个不同颜色的mesh需要重叠值范围大一点，否则mesh不能贴合到一起。
    { 12.5,22.5, 0.61, 0.81, 0.67 , 0.3 },
    { 17.5,27.5, 0.83,0.92,0.68 , 0.4 },
    { 22.5,32.5, 0.93,0.46, 0.00 , 1.0 },
    { 27.5,37.5, 0.93,0.00, 0.00 , 1.0 },
    { 32.5,42.5, 0.55,0.00, 0.00 , 1.0 },
    { 37.5,999.0,0.63,0.00, 0.63 , 1.0 },
};

int generateMeshAndTiler(
    int dim1, int dim2, int dim3,
    float k0, float b0,
    float k1, float b1,
    float k2, float b2,
    vector<float>& imgData, //size should be dim1*dim2*dim3;
    int originalMC,
    vector<MeshRange>& meshRanges,
    float reduceFraction,
    bool preSmooth,
    bool onlyLargest,
    bool fillBubbles,
    int postSmooth,
    bool verbose,
    int quality,
    const char* outTileDirWithouSlash
);


int main(int argc,char* argv[])
{
    cout << "v1.0" << endl;//version
    cout << "v1.1 add support for 126 bands grib file." << endl;
    cout << "usage:yy_era5_105_to_windmesh /some/dir/...grib /out/dir" << endl;

    if (argc != 3) {
        cout << "argc not 3" << endl;
        return 11;
    }
    string gribfile = argv[1];
    string outdirNoSlash = argv[2];

    GDALAllRegister();
    int xsize = 0;
    int ysize = 0;
    int gbands = 0;
    GDALDataset* dsin = 0;
    const int NUM_LAYER = 21;
    vector<UvLayer> uvlayers(NUM_LAYER); //数组顺序从0到大，表示海拔高度从高到地表。
    cout << "reading uv datas..." << endl;
    try
    {
        dsin = (GDALDataset*)GDALOpen(gribfile.c_str(), GA_ReadOnly);
        xsize = dsin->GetRasterXSize();
        ysize = dsin->GetRasterYSize();
        gbands = dsin->GetRasterCount();
        if (xsize == 1440 && ysize == 721 && (gbands == 105 || gbands==126) ) {
            for (int ilayer = 0; ilayer < NUM_LAYER; ++ilayer) { // 105/5 = 21
                int iu = ilayer * 5 + 4;
                int iv = ilayer * 5 + 5; 
                if (gbands == 126) {
                    iu = ilayer * 6 + 5;
                    iv = ilayer * 6 + 6;
                }
                uvlayers[ilayer]._u.resize(1440 * 721);
                uvlayers[ilayer]._v.resize(1440 * 721);
                uvlayers[ilayer]._l.resize(1440 * 721);
                dsin->GetRasterBand(iu)->RasterIO(GF_Read, 0, 0, 1440, 721, uvlayers[ilayer]._u.data(), 1440, 721, GDT_Float32, 0, 0, 0);
                dsin->GetRasterBand(iv)->RasterIO(GF_Read, 0, 0, 1440, 721, uvlayers[ilayer]._v.data(), 1440, 721, GDT_Float32, 0, 0, 0);
                for (int iy = 0; iy < 721; ++iy) {
                    for (int ix = 0; ix < 1440; ++ix) {
                        int it = iy * 1440 + ix;
                        float uu = uvlayers[ilayer]._u[it];
                        float vv = uvlayers[ilayer]._v[it];
                        uu = uu * uu;
                        vv = vv * vv;
                        uvlayers[ilayer]._l[it] = sqrtf(uu+vv);
                    }
                }
            }
        }
        else {
            stringstream ss;
            ss << "input file size not valid. xsize(1440) " << xsize << ", ysize(721) " << ysize << ", nbands(105||126) " << gbands << endl;
            throw std::runtime_error(ss.str());
        }

        if(dsin) GDALClose(dsin);
        dsin = 0;

        // 使用 550hpa - 1000hpa 层对应YYERA5的 12-20 （索引值开始为0）
        int usedLayers = 9;
        float heightPerLayer = 20000.0;
        //为了使球闭合，这里经度方向最后一列后面增加10列数据
        const int extraColumn = 3;
        if (xsize <= extraColumn + 1) {
            throw runtime_error("invalid xsize <= 4");
        }
        int usedXSize = xsize + extraColumn;
        vector<float> windLenData(usedXSize * ysize * usedLayers);
        for (int iband = 0; iband < usedLayers; ++iband) {
            float* bandPtr = windLenData.data() + iband * usedXSize * ysize;
            float* tempbufferPtr = uvlayers[20 - iband]._l.data() ;
            for (int iy = 0; iy < ysize; ++iy) {
                for (int ix = 0; ix < xsize; ++ix) {
                    *bandPtr = *tempbufferPtr;
                    bandPtr++;
                    tempbufferPtr++;
                    if (ix == xsize - 1) {
                        for (int iex = 0; iex < extraColumn; ++iex) {
                            *bandPtr = uvlayers[20 - iband]._l[iy * xsize + iex];
                            bandPtr++;
                        }
                    }
                }
            }
        }

        cout << "do meshing..." << endl;

        //气压和海拔高度换算参考
    // 1000hpa - 111 meter
    // 975hpa  - 323 meter  
    // about 212 meter per 25hpa / per layer.
    // 这里高度方向放大100倍 每层对应 ss meter 

        float isolevel = 10.0;// by qgis try.
        int isoDarkMediumBright123 = 2;
        float reduceFraction = 0.25;
        bool preSmooth = true;
        bool onlyLargest = 0;//0 for all
        bool fillBubbles = false;
        int postSmooth = 0;
        int quality = 1;
        int originalMC = 0;
        bool verbose = false;

        int ret0 = generateMeshAndTiler(
            usedXSize, ysize, usedLayers,
            // 0.25, 0.0, //from 0.0 to 360.0;
            0.25, -180.0, //from -180 to +180.0
            -0.25, 90.0,
            heightPerLayer, 0.0,
            windLenData,
            originalMC,//Original marching cubes (0=Improved Lewiner, 1=Original, default 0)
            g_MeshRanges, //大于这个值的时候认为是有效填充，反之为空气
            reduceFraction,
            preSmooth,
            onlyLargest,
            fillBubbles,
            postSmooth,
            verbose,
            quality,
            outdirNoSlash.c_str()
        );

        cout << "mesh done: " << outdirNoSlash << endl;

        
    }
    catch (exception& ex) {
        cout << "runtime error : " << ex.what() << endl;
        if(dsin) GDALClose(dsin);
        dsin = 0;
        return 12;
    }
    cout << "done" << endl;
    return 0;
}



JsonObject& createATileMeshJsonObject(DynamicJsonBuffer& jsonbuffer,
    YYTileMesh& current,
    vector<YYTileMesh>& tileArray,
    double geoError)
{
    Ecef_Tiler etiler;
    JsonObject& theObj = jsonbuffer.createObject();

    vec3d bbmin;
    vec3d bbmax;
    current._mesh.computeBbox(bbmin, bbmax);

    JsonObject& bvol = theObj.createNestedObject("boundingVolume");
    JsonArray& box = bvol.createNestedArray("box");
    box.add(bbmin.x * 0.5 + bbmax.x * 0.5);
    box.add(bbmin.y * 0.5 + bbmax.y * 0.5);
    box.add(bbmin.z * 0.5 + bbmax.z * 0.5);
    box.add(bbmax.x * 0.5 - bbmin.x * 0.5); box.add(0.0); box.add(0.0);
    box.add(0.0); box.add(bbmax.y * 0.5 - bbmin.y * 0.5); box.add(0.0);
    box.add(0.0); box.add(0.0); box.add(bbmax.z * 0.5 - bbmin.z * 0.5);
    JsonObject& contentobj = theObj.createNestedObject("content");
    char glbname[256];
    sprintf(glbname, "./tile_%d_%d_%d_%d.glb",
        current._l,
        current._x,
        current._y,
        current._z);
    contentobj["url"] = string(glbname);
    theObj["geometricError"] = geoError;
    theObj["refine"] = string("REPLACE");
    JsonArray& children = theObj.createNestedArray("children");
    for (auto it = tileArray.begin(); it != tileArray.end(); ++it) {
        //if the tile is child of current object.
        if (etiler.isChild(current, *it)) {
            JsonObject& aChild = createATileMeshJsonObject(jsonbuffer, *it, tileArray, geoError / 2.0);
            children.add(aChild);
        }
    }
    return theObj;
}



// wf 
int generateMeshAndTiler(
    int dim1, int dim2, int dim3,
    float k0, float b0,
    float k1, float b1,
    float k2, float b2,
    vector<float>& imgData, //size should be dim1*dim2*dim3;
    int originalMC,
    vector<MeshRange>& meshRanges,
    float reduceFraction,
    bool preSmooth,
    bool onlyLargest,
    bool fillBubbles,
    int postSmooth,
    bool verbose,
    int quality,
    const char* outTileDirWithouSlash
) {
    float srow_x[4] = { k0,0,0,b0 };
    float srow_y[4] = { 0,k1,0,b1 };
    float srow_z[4] = { 0,0,k2,b2 };

    short dim[3] = { dim1, dim2 , dim3 };

    //gltf
    YYMesh fullbodyMesh;
    for (int imesh = 0; imesh < meshRanges.size(); ++imesh) {
        cout << "start imesh " << imesh << endl;
        vec3d* pts = NULL;
        vec3i* tris = NULL;
        int ntri, npt;

        vector<float> copyData = imgData;
        int asize = copyData.size();
        float* dataPtr = copyData.data();
        for (int tit = 0; tit < asize; ++tit) {
            if (*dataPtr > meshRanges[imesh].valmin && *dataPtr <= meshRanges[imesh].valmax) {
                //keep the value.
            }
            else {
                *dataPtr = 0.f;
            }
            dataPtr++;
        }

        //const int isolevel = 2;// 0-nodata , 2-filled, nothing else
        if (meshify(copyData.data(), dim, originalMC, meshRanges[imesh].valmin, &tris, &pts, &ntri, &npt, preSmooth, onlyLargest, fillBubbles, verbose) != EXIT_SUCCESS)
            return EXIT_FAILURE;
        apply_sform(tris, pts, ntri, npt, srow_x, srow_y, srow_z);
        double startTime = clockMsec();

        YYPrimitive tempPrim;
        YYWriteGltf::transformLlaPrimitive2EcefPrimitive(pts, tris, npt, ntri,
            meshRanges[imesh].r,
            meshRanges[imesh].g,
            meshRanges[imesh].b,
            meshRanges[imesh].a,
            tempPrim
        );
        fullbodyMesh._prims.push_back(tempPrim);

        free(tris);
        free(pts);
    }

    //Tiler
    Ecef_Tiler eTiler;
    const double earthRadius = 6378000;//in meter
    const double maxTileRadius = earthRadius + 310000;//add  310km for tiler
    vector<YYTileMesh> tileMeshArr;
    string tileError;
    eTiler.cutTheTiles(3, maxTileRadius, fullbodyMesh, tileMeshArr, tileError);
    cout << "good tile count " << tileMeshArr.size() << endl;
    cout << "writing tiles into glb files" << endl;
    for (int itile = 0; itile < tileMeshArr.size(); ++itile) {
        char outtilename[256];
        sprintf(outtilename, "%s/tile_%d_%d_%d_%d.glb", outTileDirWithouSlash,
            tileMeshArr[itile]._l,
            tileMeshArr[itile]._x,
            tileMeshArr[itile]._y,
            tileMeshArr[itile]._z
        );
        YYWriteGltf writer3;
        writer3.writeYYMeshToGlbFile(outtilename, tileMeshArr[itile]._mesh);
    }

    //write tileset.json

    if (tileMeshArr.size() > 0)
    {
        ArduinoJson::DynamicJsonBuffer jsonbuffer;
        JsonObject& doc = jsonbuffer.createObject();
        JsonObject& assetobj = doc.createNestedObject("asset");
        assetobj["version"] = "1.0";
        doc["geometricError"] = 10000000.0;
        double geometricError0 = 200000;
        YYTileMesh& rootTileMesh = tileMeshArr[0];
        for (auto itTile = tileMeshArr.begin(); itTile != tileMeshArr.end(); ++itTile) {
            if (itTile->_l == 0) {
                rootTileMesh = *itTile;
                break;
            }
        }
        JsonObject& tileRootObj = createATileMeshJsonObject(jsonbuffer, rootTileMesh, tileMeshArr, geometricError0);
        doc["root"] = tileRootObj;

        string tilesetJsonName = string(outTileDirWithouSlash) + "/tileset.json";
        string jsontext;
        doc.printTo(jsontext);
        cout << "writing " << tilesetJsonName << endl;
        std::ofstream ofs(tilesetJsonName.c_str());
        ofs << jsontext;
        ofs.close();
    }
    cout << "done!" << endl;
    return EXIT_SUCCESS;
}

